* Program to construct bootstrap percentiles of noise-signal ratio from real-time gap
* data over the 1980-1992 sample period.
* Note:  blockwidth and nobs are set by hand, and this program requires the number
* of blocks to divide evenly into the sample for the circular MBB procedure.


* This program:  Uses a block width of 4 and the standard deviation-based noise-signal ratio


calendar 1947 1 4
allocate 2013:04+5000

seed 4892

open data table1_boot_data.rat
data(format=rats) / gb_gap_rt gb_gap_gdp94q3

* GB_GAP_RT:       Real-time gap estimate
* GB_GAP_GDP94Q3:  Final gap for 1980-1992 period


* Put the relevant data series into a vector (note that
* the make command will actually yield an array...)

make rtgap_vec 1980:01 1992:04
# gb_gap_rt
disp %rows(rtgap_vec); * Check that this equals NOBS, below

make fingap_vec 1980:01 1992:04
# gb_gap_gdp94q3
disp %rows(fingap_vec); * Check that this equals NOBS, below

compute gaprev = fingap_vec-rtgap_vec
* compute sqgaprev = gaprev.^2.0

compute blockwidth=4
compute nobs=52
compute numblks=13
compute numdraws=5000

declare vector xxx(numblks)
declare vector[integer] entries(nobs)
declare vector b_rev(nobs)
declare vector b_fin(nobs)
declare vector n_s_ratios(numdraws)



do drawiter=1,numdraws


ewise xxx(i)=%raninteger(1,nobs)

* Use these two lines of code (with NUMDRAWS=1) to test the calculations (they just return the original series)
* compute qqq=%seqa(1,blockwidth,nobs)
* ewise xxx(i)=qqq(i)


do iter=1,numblks
  do subiter=1,blockwidth
    compute entries(((iter-1)*blockwidth)+subiter)=%clock(fix(xxx(iter))+(subiter-1),nobs)
  end do subiter
end do iter

do iter=1,nobs
  compute b_rev(iter)=gaprev(entries(iter),1)
  compute b_fin(iter)=fingap_vec(entries(iter),1)
end do iter

compute sqdev_r=(b_rev-%avg(b_rev)).^2.0
compute nnn_r=%rows(sqdev_r)
compute nscale_r=float(nnn_r)/float(nnn_r-1)
compute stddev_r=sqrt(nscale_r*%avg(sqdev_r))

compute sqdev_f=(b_fin-%avg(b_fin)).^2.0
compute nnn_f=%rows(sqdev_f)
compute nscale_f=float(nnn_f)/float(nnn_f-1)
compute stddev_f=sqrt(nscale_f*%avg(sqdev_f))

compute n_s_ratios(drawiter)=stddev_r/stddev_f

end do drawiter



set sd_nsr_series 1 numdraws = n_s_ratios(t)
statistics(fract) sd_nsr_series

